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A Closed-Form Formulation of HRBF-Based 

Surface Reconstruction 

Shengjun Liu, Charlie C. L. Wang, Senior Member, IEEE, Guido Brunnett, and Jun Wang 


Abstract —The Hermite radial basis functions (HRBFs) implicits have been used to reconstruct surfaces from scattered 
Hermite data points. In this work, we propose a closed-form formulation to construct HRBF-based implicits by a quasi-solution 
approximating the exact solution. A scheme is developed to automatically adjust the support sizes of basis functions to hold the 
error bound of a quasi-solution. Our method can generate an implicit function from positions and normals of scattered points 
without taking any global operation. Working together with an adaptive sampling algorithm, the HRBF-based implicits can also 
reconstruct surfaces from point clouds with non-uniformity and noises. Robust and efficient reconstruction has been observed in 
our experimental tests on real data captured from a variety of scenes. 

Index Terms —Hermite Radial Basis Functions, Quasi-solution, Closed-Form, Surface Reconstruction 
- ♦ - 


1 Introduction 

R ECONSTRUCTING surface from a set of unorga¬ 
nized points equipped with normal vectors is 
an important topic in various fields such as com¬ 
puter graphics, reverse engineering, image process¬ 
ing, mathematics, robotics and CAD/CAM. A lot of 
research approaches have been devoted to develop 
surface reconstruction methods, in which implicit sur¬ 
face fitting based on Radial Basis Functions (RBF) is 
successful in dealing with noisy and incomplete data 
(e.g., [l]-[3]). 

Recently, implicits based on Hermite Radial Basis 
Functions (HRBF) were presented to interpolate data 
points to the first order [4], It is robust and effec¬ 
tive to deal with coarse and non-uniformly sampled 
points, close surface sheets, and surfaces with fine 
details. However, interpolating both positions and 
normals of points leads to the computation of solving 
a 4n x 4 n linear system for an input with n points. It 
is impractical due to the expensive computation. The 
system becomes sparse when the Compactly Supported 
Radial Basis Functions (CSRBF) are used as the kernel 
functions. However, this attempt on improving the 
efficiency can bring in a more challenging problem 
- numerical stability. 

HRBF-based quasi-interpolation is presented in this 
paper to overcome the computational problem. Quasi- 
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Fig. 1. The method proposed in this paper can effi¬ 
ciently reconstruct a surface from a set of noisy and 
incomplete points - e.g., the indoor scene shown here 
with 922k points. Our reconstruction takes only 60 sec¬ 
onds to generate a mesh surface that has the similar 
quality as the state-of-the-art [5] but the computation is 
7.6x faster. The parameter, s = 3.5, is employed in this 
reconstruction. 


interpolation is a kind of approximate interpolation 
that fits implicits by weighted averages of the values 
at given points. The most attractive property of quasi¬ 
interpolation is to reconstruct a surface from a set of 
points without solving linear systems - i.e., with a 
closed-form formulation. This can make the computa¬ 
tion of HRBF-based surface reconstruction stable and 
efficient. As shown in Fig.l, the mesh surface can be 
efficiently reconstructed from an input set with 922fc 
points by our method in 60 seconds. Comparing to the 
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recently published Floating Scale Surface Reconstruction 
(FSSR) that also avoids applying global operations, 
our method is about 7.6 x faster, s = 3.5 is the only 
parameter used in our approach. Here, s is defined as 
an amplifier of CSRBF kernels' support sizes, which 
control the maximal number of data points covered 
by each support (details can be found in Section 4.1). 
Moreover, we have analyzed the error-bound between 
our closed-form solution and the solution obtained by 
solving linear systems (see Section 3.2). Specifically, 
the error-bound exists when the number of points 
covered by the support of each kernel is capped by 
a fixed number. To overcome the problem caused 
by high non-uniformity on points, an algorithm is 
developed to select a sub-set of given points to serve 
as centers of kernel functions, which is an optional 
step in our framework of reconstruction. 

1.1 Main result 

In this paper, we propose a closed-form formulation 
for computing the quasi-solution of HRBF-based sur¬ 
face reconstruction from scattered data points. 

• The method can construct a signed scalar func¬ 
tion by directly blending the positions and nor¬ 
mals of points without any global operation. The 
computation based on CSRBF is local and robust. 

• Errors between the quasi-solution and the exact 
one are bounded after applying an automatical 
scheme to adjust the support sizes of basis func¬ 
tions. 

• Combining with an adaptive center selection 
algorithm, surface reconstruction based on our 
method can remove the artifacts resulted from the 
non-uniformity and noises. 

As a local approach, our method is efficient and 
scalable. This is well-suited for highly parallel imple¬ 
mentation as well as distributed/progressive recon¬ 
struction. Note that, the compactly-supported basis 
functions results in open meshes and leaves holes in 
the region does not have enough number of points, 
which fits the application of reconstructing outdoor 
scenes very well. 

1.2 Related work 

The problem of surface reconstruction from point 
cloud has been studied in literature for more than 
two decades. A comprehensive review of all these 
works has beyond the scope of our paper. More 
discussion and comparison on different surface recon¬ 
struction methods can be found in [6]. We only give 
an overview of implicit function based reconstruction 
with a focus on RBF-based formulations. 

After using signed distance field in [7] to recon¬ 
struct mesh surface from point clouds, implicit func¬ 
tions have gained a lot of attention in surface re¬ 
construction because of its ability to handle topo¬ 
logical change and fill holes. Example approaches 


include RBF-based methods [1], [2], [4], [8]—[16], Pois¬ 
son surface reconstructions [17], [18], smooth signed 
distance method [19], moving least-squares based 
methods [20]—[25], wavelets based method [26], and 
Partition-Of-Unity (POU) based methods [27], [28]. The 
methods based on RBF implicits are popular for their 
capability of handling sparse point clouds. Generally, 
RBF-based methods transform the reconstruction into 
a problem of multi-variational optimization, where 
enforcing the interpolation constraints results in a 
linear system. Solving the linear system is an im¬ 
portant but time-consuming step for the RBF-based 
reconstruction. To obtain a non-trivial solution, RBF- 
based methods usually require the provision of extra 
offset-points (ref. [1], [2]) that can be obtained by 
shifting data points along their normal directions. 
However, it is not easy to find an optimal value for 
the offset scalar. The positions of these offset points 
is also difficult to determine, especially when the 
scanned model has thin sheets and the distribution of 
input points is irregular. To avoid generating offset- 
points, Ohtake et al. [9], [13] used a signed function 
which includes basic approximations and local details. 
The basic approximation formed by local quadratic 
functions with POU is considered as an offset function 
which constructs the non-zero constraints for fitting 
local details with RBFs. A very important information 
that describes shape to the first order - the normal 
vectors of a model have not been well utilized in these 
approaches of surface reconstruction. 

Some prior works (e.g., [10], [15]) deduced from 
the statistical-learning perspective avoid generating 
offset-points in surface reconstruction, where normals 
were directly used in a variational formulation. Re¬ 
cently, Macedo et al. [4] derive an implicit function 
from the Hermite-Birkhoff interpolation with RBFs. 
They enhance the flexibility of HRBF reconstruction 
by ensuring well-posedness of an approach combin¬ 
ing regularization. However, given a set with n points 
and n normal vectors, these methods give a 4 n x 4n 
linear system to be solved, which limits the number 
of points can be involved in the reconstruction. 

Quasi-interpolation is a method in the field of 
function approximation. It is simple, efficient, and 
computational stable. In the early work of quasi¬ 
interpolation [29], a function approximating a given 
data set is defined by a weighted average of the values 
at the data points. The idea has been used in [30] for 
surface reconstruction. The quasi-interpolation with 
radial basis functions has been studied in [31]. Re¬ 
cently, Wu and Xiong [32] developed a new method 
to construct kernels in a quasi-interpolation scheme 
by the linear combination of scales, where the kernels 
are in the form of RBF. Han and Hou discussed 
quasi-interpolation by RBF and suggested values of 
the shape parameters in [33], which is a constructive 
method for obtaining a family of quasi-interpolations. 
Liu and Wang generalized the regularly sampled grid 
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points to 3D scattered points in [34], They proposed a 
multi-level quasi-interpolation method which is based 
on POU and the RBF method. However, linear sys¬ 
tems still need to be solved in their method. Locally 
supported basis functions satisfying the property of 
POU are used in the recent work FSSR [5], in which 
weighted average is employed to fit implicit functions 
to the input set of points. However, there is no error- 
bound guaranteed in their approach. Experimental 
tests show that our method can generate results with 
similar quality as FSSR but has 7.61 to 98.3 times 
speedup. 

There are schemes for finding a subset of 'optimal' 
centers from the point set in order to obtain fast recon¬ 
structions. Carr et al. [1] proposed a greedy algorithm 
that iteratively appends centers which are correspond¬ 
ing to the maximal residual of the current RBF fitting 
until a desired accuracy is reached. Samozino et al. 
presented the reconstruction with voronoi centered 
RBFs [11]. Ohtake et al. [13] proposed a reconstruc¬ 
tion method which combines an adaptive POU ap¬ 
proximation with least-squares RBF fitting. Different 
from [13], we adopt a quadric error function based 
on positions and normals of a point and its neighbors 
instead of local quadratic approximation. Therefore, 
our center selection step can generate a good spherical 
cover in a non-iterative way. The selected centers of 
spheres describe the input shape with a small error. 

The rest of our paper is organized as follows. We 
first introduce the surface reconstruction based on 
regularized HRBF implicits in Section 2. Section 3 pro¬ 
vides our formulation in the closed-form and derives 
the error bound of our formulation. Section 4 presents 
the algorithms for different steps of reconstruction, in¬ 
cluding parameters tuning, isosurface extraction and 
center selection (optional). After that, the results of 
experimental tests are shown and discussed in Section 
5. Lastly, our paper ends with the conclusion. 

2 HRBF Implicits 

The HRBF implicits [4] are built upon the theory 
of Hermite-Birkhoff interpolation with radial basis 
functions [3]. In this section, we briefly describe how 
to use HRBF implicits to solve the problem of surface 
reconstruction from scattered points. 

Definition 1 Given a set of data V = {pi, P 2 , • ■ ■ , p„} 
with unit normals AT = {n 1 ,n 2 ,--- . n„}, the HRBF 
implicits give a function / interpolating both the 
points and the normal vectors as 

n 

/(x) = WM* ~ Pi) - (bj, v<p(x - Pj)>}, (1) 

i =1 

where <p : SR 3 1 —i 5ft is defined by a radial basis function 
tp(x) = </> p (||x||), (•,•) denotes the dot-product of two 
vectors, and V is the gradient operator. 


The scalar coefficients, aj € 3?, and the vector coeffi¬ 
cients, b ; G -ft 3 , can be determined by the constraints 
of interpolation as 

f(Pi) = c and V/ (p,;) = n,, (z = 1, 2, • • • , n) (2) 

with c being a constant value for the implicit function, 
c = 0 is used for surface reconstruction. Applying 
the constraints (2) to Eq.(l), we obtain a linear system 
with equations 

E"=i WMPi - Pi) - (bj, Vy>(pi - Pj))} = c, 

E"= 1 {ajVip{pi - p j) - bjHpfp, - Pj )} = n i5 

where i = 1, 2, • • • ,n and H is the Hessian operator 
applied on <p(-). The linear system can be rewritten in 
a matrix form as 

AA = y, (4) 

where A and y are 4n vectors with the z-th blocks 
being [di,b;] T and \c. n,] r respectively. Here, A is a 
4 n x 4n coefficients matrix which are assembled from 
n x n blocks. Each block A ij is a 4 x 4 sub-matrix 
corresponding to a pair of RBF centers (p ? , p ; ). 

-A- — (A ij'jn.xm 

A = (<p(Pi-Pj) -V<p(pi - PjA (5) 

" P ->) _I M P < ~ p J )/ 4 x 4 " 

In this paper, we use a Wendaland's CSRBF [35] as 
the kernel function 

<M0 = 4>{r/p) 

(l-f) 4 (4f + l), t e [0,1], (6) 

0, otherwise, 

where p is the support size, and r is the Euclidean 
distance between a query point and the center of a 
RBF. Note that, different support sizes can be used at 
different centers. Solving Eq.(4), an implicit function 
/(x) can be determined at the space spanned by the 
supports of centers {p,}. To make the system matrix of 
RBF interpolation better conditioned, a regularization 
term with coefficient 17 is usually added when using 
RBF interpolants to solve the surface reconstruction 
problem (ref. [36]). That is, 

(A + z?I)A = y. (7) 

An example is shown in Fig.2 to demonstrate the 
effectiveness of regularization. 

3 Formulation 

This section provides a closed-form formulation for 
solving the HRBF-based surface reconstruction prob¬ 
lem via quasi-solution. Error-bound of the approxi¬ 
mation is also derived. 
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Fig. 2. Surface reconstructions by the HRBF implicits 
with (left) and without (right) the regularization term. 
Artifacts will be produced when the regularization term 
is not added. In this example, 77 = 100/p 2 is employed. 


3.1 Quasi-solution in closed-form 

When increasing the number of centers in the linear 
system of FtRBF-based surface reconstruction (i.e., 
Eq.(7)), numerical instability and expensive compu¬ 
tation become intensively remarkable. Here, we in¬ 
vestigate a closed-form formulation derived from the 
theorem of quasi-interpolation to reconstruct surface 
in a more stable and efficient way. 

Quasi-interpolation technique can reconstruct a 
function interpolating a given data set by computing 
weighted averages of the values at the data points 
[29]. Specifically, considering an exact interpolant 
s(x) = with the constraints gfc) = f , 

of function values, the function g(x) can be well 
approximated by letting A i = /,. That is a quasi¬ 
interpolation, g(x) = ./idifx). However, the quasi¬ 
interpolation technique cannot be directly applied 
here as our interpolation constraints consist of both 
the values and the gradients of functions (see Eq.(2)). 

Basically, we need a closed-form formulation to 
approximate the solution of Eq.(7). By means of the 
matrix computation, the quasi-interpolant with \i = 
/, is actually a quasi-solution when the coefficient 
matrix is approximated by an identity matrix I. Here, 
a similar approximation is employed in the HRBF- 
based reconstruction problem. For a CSRBF 
when there is no other center falling into the space 
spanned by its support pt, the coefficient matrix is 
degenerated from i of Eq.(5) into 

^ , . 20 20 20 , T ^ f / \ 

Du = diag( 1, -^,^,^) + 77 I 4 , Djj = 0 (t A j). 

Pi Pi Pi 

( 8 ) 

If the scenario of not containing other centers hap¬ 
pens at all CSRBF kernels, the linear system to be 
solved in Eq.(7) is degenerated into DA = y with 
D = (D ij) nxn - This leads to an approximate solution 


of Eq.(7) as 

A =T>~ l y 

_ r c pjn 1 
f 1+7 }’ 20 +7}p\ ’ ‘ ■ 


c pjnn ~j 

1-1-77 ’ 20+77 p'i t 


(9) 


The zero level-set is usually employed in surface 
reconstruction (i.e., c = 0 in Eqs.(2) and (3)). c = 0 
is used in all formulas in the rest of this paper. As a 
result, the coefficients of the i-th basis function can be 
approximated by 


p? 

ai ~ 0 and b, « ——^^ 
20 + ppf 


which give an approximate function of /(x) in a 
closed form as 


/« = -£< 

j =1 


Pj 

20 + 77 p: 


rnj,Vp(x-p j)). 


( 10 ) 


By this implicit function, we can apply the polygo- 
nization techniques (e.g., [37], [38]) to tessellate the 
isosurface of /(x) = 0 into a polygonal mesh - 
the surface is reconstructed from scattered Hermite 
points. 


3.2 Error-bound analysis 

The error between the quasi-solution A and the exact 
solution A of Eq.(7) must be bounded to make the 
closed-form formulation useful. The analysis is given 
below. 


Lemma 1 Defining AA = (A + 77 I) — D and AA = 
A — A, the error of approximation is bounded as 


||AA||oo < 


D 


-n 


IA All 


when IID 


1 — ||D _ 1 ||c 5 o II AA|| 
IIAAIU < 1 . 


D 


- 1 , 


( 11 ) 


Proof: By Eq.(7), we have 


(D +AA)(A +AA) = y. 

With the quasi-solution that A = D _1 y, this equation 
can be converted to 


AA = D^ 1 [—(AA)A - (AA)(AA)]. 

Then, we apply a maximum norm and get 

IIAAIU < l|D- 1 ||oo(l|AA|| 00 ||A || 00 + IIAAIU IIAAIU), 

which is also 


(1—||d- 1 |uI|aa|U)I|aa || 00 < I|d _ 1 IIoo||aa|ui|a|U- 

By the given condition ||D~ 1 || 00 ||AA|| 00 < 1, we have 
IIAAH < ll D 1 llooIIAA|U II Q| 

II A H- - ! _ || D -1|UI| AA |U 11 °°' 

Combining with Eq.(9), the lemma has been proved. 
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Assuming there are at most m other centers falling 
in the support region for each kernel, the error-bound 
of our quasi-solution can be achieved on the Wends- 
land's CSRBFs (i.e., Eq.( 6 )). 


Lemma 2 When Wendsland's CSRBFs are used, if 
their support sizes p t <E [p min , Pmax] (with p max < V 20) 
and each support region contains at most m centers 
of other CSRBFs, the error of HAAH^ is bounded by 
a constant when 


V > 



5 35 

- 1 - 2 - 

4p m in Pmin 



( 12 ) 


Proof: By the definition of the diagonal matrix D in 
Eq.( 8 ), 


ll D_1 || 


max { 


1 + ?7 ’ 20 + r/pj 


}• 


The upper bound of j|D 1 y || 00 can also be obtained 
from Eq.(9) as 


IID-Vll 


IIAIloo 


max {0, 


p 2 n x p 2 n v - 

r’j J ^3 3 

20 + rjp 2 ’ 20 + rjp 2 ’ 


P 2 3 n 3 
20 + ?y p: 




TABLE 1 

Derivatives of Wendsland’s CSRBF [35] 



f Here, x = (. x, y, z), r = y/(x - Xi) 2 + {y - y t ) 2 + (z - Zi) 2 and 

pi is the support size of the radial basis function (p;( x )- 


Here, superscripts denote the x-, y- and .-'-components 
of a vector in 3? 3 . When p t > pj > 0 and 77 > 0, 

pi = 1 > p] = 1 

20 + r/p 2 20/ p 2 + rj 20 + 77 p 2 20/ p 2 + 77 ’ 

As a result 


||D“ 


< max( 


p 2 

r max 


1 + 77 ’ 20 + 77 


|D Vlloo < 


20 + r/p. 


,2 

max 


When pj < 77 lnax < V20, we can further obtain 


ll D_1 || 


1 

I + 77' 


Now we derive the upper bound of jjAAHoo. From 
Eqs.(5) and (8), using <pjj to denote <p,;(pj) = ip(p :l — 
Pi), we can also have 


aa 

|| 00 — 

max { 

j=l,...,n 




m rx 

X^I^Cil + l g x 1 + 

i 


d+i 

dz 

^ 1 ), 

m 

Ed 

i 

d Pip 
dx 


dxdy 

l + l 

d 2 Vi,j 

dxdz 

m 

Ed 

i 

d<Pi t j 

Oy 

1 + 1 dxdy 

, 1 d2( PiJ 
' dy 2 

l + l 

d 2 <Pij 

dydz 

m 

Ed 

i 

dz 


, 1 d 2 Vi,j 

1 dydz 

l + l 

d 2 ipi,j 

dz 2 


I)}- 


By the derivatives listed in Table 1 and their corre¬ 
sponding upper bounds listed in Table 2, we can have 

1 5 5 35 

IIAAHoo < max{?77.(1 + -—),m(--h -^)} 

3 4 Pj 4 Pj p 2 


15 5 

< max{m(l + - -),m(-- 

4p m in j 


45 -,,. 


When p min < p n 
simplified to 


< a/ 20 / it can easily be further 


||AAHoo < 777( 


4p n 


35 

Pmin 


) = A. 


Summarizing all the analysis together, we have 

Apl 


IIAAIIoo < 


9 2 

J max 


(I + 77-A)(20 + 7?p^ ax ) 


(13) 


when ||D _1 || 00 ||AAjloo < A/{ 1 + 77 ) < 1. To hold this, 
it should have 1 + 77 > A — that is Eq.(12). The lemma 
has been proved. 


Remark. The requirements of, 

1) all the CSRBFs have their support sizes within 
the interval [p min , Pmax], and 

2 ) there are at most 777 centers falling in the support 
of any others CSRBF, 

can be achieved by a carefully designed parameter 
tuning algorithm (see Section 4.1). After determining 
the support sizes, the value of 77 can be chosen to hold 
Eq.(12). By scaling all models into a bounding box of 
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TABLE 2 

Error Bounds of Derivatives 


IviMI 

< 1 

|9^i(x) | 9<pi{x) | d<pi(x) | 

dx dx dx 

< — (1- -) 3 - 

Pi Pi Pi 

5 

< - (with r = Pi-) 

4 Pi 

l'TVi( x ), ,d 2 ipi{yi) 3 2 <ft(x) 

ox z oy z oz z 

20 . r n9 . r . 

< -^(1 — ) 2 (i + 2 —) 

pi pi pi 

20 

< — (with r = 0) 

Pi 

, 9 2 <pi(x) a 2 ^i( x ), , d 2 ipi{x.) 

1 *-v ex Ml 0 Ml ex M 

oxoy oxoz oxoy 

pi pi pi 

< il ( W ith r = f) 

_ lA _ 


tThe analysis is based on \x — aq| < r, \y — yi\ < r and 
1 2 — Zi\ < r, and the bound is derived by using the inequality of 
arithmetic and geometric means. 


[—1, l ] 3 € SR 3 , the algorithm in Section 4.1 can also 
enforce p max < y/20. Numerical errors generated in 
our experimental tests are discussed in Section 5.3, 
which further verify the error-bound of our method. 

4 Reconstruction Algorithm 

Given points and their normals in the input sets V and 
Af, the implicit function /(x) defined in Eq.(10) can be 
evaluated in the supported regions as an approxima¬ 
tion of the HRBF reconstruction. Specifically, with the 
help of tessellation techniques, zero isosurface of the 
implicit function /(x) can be converted into a mesh 
surface. A scheme is also developed to determine the 
support sizes of HRBFs and the coefficient r) according 
to Lemma 2 to guarantee the existence of error-bound. 

Our reconstruction algorithm consists of a scheme 
to tune parameters according to the analysis for error- 
bound and an efficient method akin to DC to extract 
the zero isosurface as a polygonal mesh in the sup¬ 
ported regions of /(x). When processing highly non- 
uniform point sets, an optional step of center selection 
is needed to reduce the artifacts caused by the non¬ 
uniformity. After that, the input sets V and A F are 
reduced into a set C = {ci,ni,--- ,c;,n?} with less 
Hermite points (i.e., I < n). 

4.1 Parameters tuning 

The support sizes of HRBF implicits, pjS, and the 
coefficient of regularization, //, should satisfy the con¬ 
dition in Eq.(12) for the existence of error-bound. 
Specifically, the following factors must be considered: 

• Pmax = niax{/9j} < \/20 gives the upper bounds 
of py 

• According to Eq.(12), with larger p rn ; n , users will 
have more flexibility to choose the value of r\ (i.e.. 


has smaller value for the right of the inequality 
in Eq.(12)). 

• On the other aspect, if increasing the support size 
Pj of a RBF centered at p j, more other centers 
will be covered by the supporting region. This 
increases the value of m in Lemma 2 and thus the 
right-hand side of the inequality in Eq.(12). 

Based on these reasons, we develop a scheme below 
to determine the values of pjS and 77 . 

Each support must cover enough number of data 
points to generate a span of the local shape; mean¬ 
while, it cannot be too large. For this purpose, we 
first construct an octree to split the input points into 
different nodes by keeping similar number of points 
in each leaf-nodes, d is then set as 3/4 of the average 
diagonal length of the leaf-nodes. The support sizes 
are temporally set as sd with s being an amplifier 
- s = 1.0 is employed for clean data sets. We then 
count the number of points covered by the supporting 
region of each CSRBF, ^> Pj . (||x — p 7 ||). The value of 
to is selected as the maximal number of data points 
covered by each of these temporary supports. After 
determining the value of m, pj of each CSRBF is en¬ 
larged by an incremental procedure until the support 
contains more than m points. By the fixed support 
sizes, the value of 77 can be chosen by Eq.(12). In all 
our tests, we assign it with a value slightly larger (e.g., 
10--) than the right-hand side of Eq.(12). 

The values of p 7 s determined by the above method 
work well on clean data but may fail in highly noisy 
input. Further tuning is needed. First of all, we enlarge 
the temporal support sizes by using s > 1.0 to 
enhance the effectiveness of denoising. This results in 
a larger to. More geometric details can be preserved 
with a smaller support size while a larger support size 
leads to a smoother reconstruction. Moreover, among 
all the PjS obtained by the aforementioned method, 
the minimal value of pjS is selected as the final sup¬ 
port size for all CSRBF kernels - i.e., uniform support 
size is adopted for highly noisy input. Examples using 
different amplifiers can be found in Fig.3. 

4.2 Efficient isosurface extraction 

Our surface reconstruction method only evaluates the 
function values of the implicit function /(x) during 
the isosurface extraction step. An variation of DC 
algorithm [38] is developed to extract zero isosurfaces 
from the regions spanned by the RBFs in Eq.(10). In 
our method, the reconstruction with compactly sup¬ 
ported implicit functions leads to open mesh surfaces 
and leaves holes in regions not covered by the sup¬ 
ports of RBFs. This is very useful for reconstructing 
the scenes that have not been completely captured 
(e.g., the scene shown in Fig.l). Moreover, as /(x) 
is defined in a closed-form, the function evaluation 
(therefore also the isosurface extraction) is highly 
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Fig. 3. Results with different regularization can be obtained by using different amplifiers on a noisy point set: 
(from left to right) input noisy points and our reconstruction results. In these examples, 77 is chosen as a value 
slightly greater than the right of Eq.(12). 



Fig. 4. Adaptive HRBF implicits are generated by our method with the help of center selection: (a) the input set 
with 100,371 points in high non-uniformity, (b) the reconstruction using all points as centers of HRBF implicits will 
easily lead to holes in the sparse regions, (c) the spherical cover - the spheres are displayed in radii as 1/4 of the 
real ones, (d) the selected 13,446 centers of RBFs, and (e) reconstruction from the selected centers - no hole 
will be generated as the densities of centers in the left and the right are similar to each other. The support sizes 
are determined by the method in Section 4.1 with (b) s = 2.1 and (e) s = 1.0, both of which lead to m = 120. 


scalable and can be efficiently performed by local 
operations. 

Voxels with a fixed width w are constructed and 
those intersecting the isosurface /(x) = 0 are first 
searched in the supported regions around the centers 
{pj}. The voxels are constructed only when 1) all of its 
eight corners have function values defined and 2) the 
function values have different signs. For each edge e of 
a voxel, the intersection between e and the isosurface 
can be determined by a bi-sectional search when two 
endpoints of e have different signs in /(•). The normal 
of this intersection point is assigned as V/. For each 
voxel, a vertex is constructed at the position v that 
minimize the quadratic-function ff , f( v ~ q^) • n, i( ) 2 , 
where q, and n q . are the intersection points on edges 
and the normal vectors at the intersections. For each 
edge e with intersection, a quadrangle is constructed 
by linking the vertices in the four voxels around e. As 
a result, the final surface mesh can be obtained. 

In our current implementation, the mesh surface 
is extracted from voxels with a fixed width. More 
sophisticated algorithm can be developed to construct 
an octree based hierarchy to extract triangles adap¬ 
tively on a partial region of implicit surfaces (e.g., 
[39]). All the operations in our isosurface evaluation 


and extraction steps are local. It is easy to be imple¬ 
mented in parallel on many-cores or in a distributed 
environment. By the algorithm's locality, a progres¬ 
sive reconstruction can be easily developed by only 
updating a local region when new points (as centers 
of CSRBFs) are added. 

4.3 Center selection 

This is an optimal step to be applied when high 
non-uniformity is observed on the input points. For 
such a point set, the direct reconstruction by using all 
points as centers of CSRBFs results in a reconstruction 
with holes in the sparse regions (see Fig.4(b) for an 
example). Flere, we adaptively select samples from V 
and AT to form a subset C. The Hermite points in C 
will be used as centers of FIRBF in the above method 
to obtain a better surface reconstruction. Each center, 
Cfc, is also associated with a radius, rp., which therefore 
forms a local spherical cover of the given points. 

Definition 2 The degree of coverage (DoC) at a point 
x £ SR 3 is defined as a function 

1 

s( x ) = Cfe ID 

k=1 


( 14 ) 
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according to a set of down-sampled Hermite points, change. As a result, a better reconstruction can be 
C = {ci, ni, ■ • • , c;, u/}. obtained (see the result shown in Fig.4(e)). 


We wish to generate a minimal spherical cover by 
controlling DoC in an iterative procedure. 

The basic idea of center selection is to form spheri¬ 
cal covers by letting DoC at every point in V not less 
than a criterion g min (i.e., Vpj G V, g(pj) > g min ). To 
this end, the following steps are iteratively run until 
the criterion is satisfied at all points. 

1) In the initial step, C = 0 and <jj = 0 is assigned 
to all points p , G V. 

2) Randomly selecting vj points with their DoCs less 
than g m i n . Among these w points, the point with 
the smallest g(-) is chosen as a center c/ i: to add 
into C together with its normal vector. 

3) The radius r k of sphere centered at c k is then 
determined by a quadric-error function 


q(c k ,r k ) 


Ej - Cfc||)(nj • (c fc - p j)) 2 

Ej^^dlPj -cfe||) 


which evaluates how curved the surface inside 
the sphere is - the shape is represented by sample 
points in V. In other words, for a highly curved 
region, a sphere with smaller r k should be used 
to reduce the error. Here, 8 j is the average of 
the squared distances between a point p j and 
its 15 nearest neighboring points. The value of 
8 j indicates a weight of point density. The bi¬ 
sectional search is taken to obtain a maximal r k 
that still satisfies 

q(c kl r k ) < q e rr L 

with L being the diagonal length of the input 
points' bounding box. 

4) Updating DoC at all points p j within the range 
11Pj - Cfc|| < r k while g(pj) < g min . DoC of c k 
is assigned as g m i n to avoid being selected as 
candidates of centers once again. 

5) Go back to step 2) until DoC at all points are not 
less than g min . 

Note that, this iterative procedure is a variant of our 
prior work in [40] with certain modification to fit 
the formulation of CSRBF. The efficiency of compu¬ 
tation has also been improved. An example result of 
our minimal spherical covering is shown in Fig.4(c), 
where the selected centers of RBFs are displayed as 
spheres. Colors are used to represent the sizes of 
spheres with red for the smallest and blue for the 
biggest ones. Samples adaptive to the geometric de¬ 
tails have been illustrated as Fig.4(d). In all examples 
of this paper, g m ; n = 1.5, q err = 5 x 10” 4 and w = 15 
work well. 

With the centers selected above, a reduced implicit 
function can be obtained by Eq.(10) but with fewer 
centers of CSRBFs. After applying the center selection 
step, the density of centers at each region becomes 
compatible to its neighboring regions - i.e., no sharp 


5 Results and Discussion 

An surface reconstruction algorithm based on the 
closed-form formulation of HRBF implicits has 
been implemented with Microsoft Visual C++ and 
OpenGL. We evaluate our methods on a PC with two 
Intel Core i7-2600K CPUs at 3.4GHz plus 16GB RAM. 
Our approach has been applied to various data with 
up to fourteen millions of points (the surface can be 
reconstructed in 78.9 sec.). Our results are compared 
with a variety of approaches - see the examples and 
discussions presented below. All the models are re¬ 
scaled into a bounding-box of [—1, l] 3 G 9i 3 . 

5.1 Comparisons 

Firstly, we test the performance of our approach on 
sets of clean points, which are uniformly sampled 
from polygonal meshes. Four models, Ramesses, Rap¬ 
tor, Momento and Neptune, are sampled into sets 
with 0.58M ~ 4.98M points. Our results are compared 
with three prior methods, including the Multiple Par¬ 
tition of Unity (MPU) reconstruction [28], the Smooth 
Signed Distance (SSD) reconstruction [19] and the 
Screened-Poisson reconstruction [18]. Comparisons are 
shown in Fig.5. We employ the 10-th depth of octree 
in the SSD and Screen-Poisson methods to generate 
results in Fig.5. For MPU and our method, we adjust 
the resolutions of polygonization methods to extract 
meshes with similar numbers of triangles as SSD and 
Screened-Poisson. The parameter Max Error of MPU 
is set as 0.001 times of the model's size. Default values 
are used for other parameters. From observation, it 
can be found that geometric details on the original 
mesh can be well preserved by our method while 
being smoothed out in some prior methods. Pub¬ 
licly available software, Metro tool [41], is employed 
to compute the average shape approximation error 
between the reconstructed surface and the original 
mesh. A bar chat of errors is given in the upper- 
right of Fig.5. Our method can always generate more 
accurate results than SSD and MPU. Meanwhile, our 
results have similar accuracy comparing to Screened- 
Poisson. 

Table 3 gives the computational statistics of tests 
on these models. Due to the close-form formula¬ 
tion, our method does not need any global operation 
such as solving a large linear system. Therefore, its 
computational time is only spent on constructing an 
octree to computing the support size and the step 
of function value evaluation in iso-surface extraction. 
Both SSD and Screened-Poisson need to solve linear 
systems globally. In Poisson reconstruction, the multi¬ 
grid solver performs a constant number of conjugate- 
gradient iterations at each level, which gives linear 
complexity w.r.t to the number of nodes in the octree. 
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Fig. 5. Experimental tests on clean point cloud that is uniformly sampled from fours mesh models - (a) Harnesses 
(0.580 M pts.), (b) Raptor (l.OOM pts.), (c) Momento (2.52 M pts.) and (d) Neptune (4.98M pts.). For illustration, 
only 1/10 points are displayed for the points of the first three models, and 1/20 points for the Neptune model. 
The reconstructions by different methods including SSD, MPU, Screened-Poisson and ours are shown and 
compared in the left. A bar chat is also given to report the average shape approximation errors on different 
reconstructions by the Metro tool [41], We use sd = 0.003 as the support size and w = 0.001 for the grid width 
of polygonization in all examples here. To conduct a fair comparison, similar number of triangles are generated 


through the polygonization for different approaches. 

The SSD reconstruction uses conjugate-gradients to 
determine all the coefficients simultaneously, which 
has a complexity of 0(n 15 ). This leads to a sig¬ 
nificantly slower performance on models with large 
number of points (see Table 3). In MPU reconstruction, 
only local fitting is taken at leaf-nodes of an octree. 
These surfaces are blended together to form the resul¬ 
tant surface, which is fast but still slower than ours. 
Moreover, our method generates results with smaller 
shape approximation error than MPU (see Fig.5). In 
summary, our method is the fastest method and can 
generate similar results as the best of other three in 
terms of quality. 


TABLE 3 

Runtime performance of different reconstruction 
approaches on clean point sets! 


Model 

Pts. 

Time in 

Seconds* 


SSD 

MPU 

Poisson 

Ours 

Ramesses 

0.58M 

14,314 

61.2 

40.8 

8.3 

Raptor 

l.OOM 

1,799 

47.2 

31.6 

6.8 

Memento 

2.52M 

24,195 

138.8 

92.6 

20.4 

Neptune 

4.98M 

6,772 

139.4 

114.0 

18.9 


‘Note that, the time reported here includes both the 
surface reconstruction and the mesh extraction. 

Mo have a fair comparison, similar number of triangles 
are generated for different approaches. 


5.2 Raw data 

In practice, the real data obtained from an acquisi¬ 
tion process are usually large and noisy point sets. 
Meanwhile, the data sets are incomplete in most cases 
(e.g., the Aquarius and the Horse models shown in 
Fig.6). The recently developed FSSR method [5] is 
targeting on fast surface reconstruction from such 
kind of real data. We compare our method with FSSR 


on two sets of real scanned data in Fig.6. Fuhrmann 
and Goesele [5] assumed the scale of an input point 
set is known, which however is not the case here. 
Although it can be computed by using the average 
distance to /{'-nearest neighbors, it is difficult to set a 
proper k for reconstructing a smooth surface. In order 
to make an appropriate comparison, we use 1/3 of 
the average support size determined by our method 
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FSSR 

Ours 


Num. of Points 

Num. of 

Time in Seconds 

Num. of 

Time in Seconds 

Model 

After Consolidation 

Triangles 

One-core 

8-cores 

Triangles 

One-core 

8-cores 

Aquarius 

833.4k 

624,916 

826.9 

172.3 

624,363 

82.5 

19.5 

Horse 

239.8k 

241,614 

263.7 

56.4 

623,772 

37.3 

9.8 


Fig. 6. Examples of surface reconstruction on incomplete set of points: (a) reconstruction from raw data and (b) 
reconstruction from data sets processed by the consolidation method [42], Our results are comparable with that 
obtained by FSSR but ours is 5.76x ~ 10.Ox faster. 



Input Screened-Poisson FSSR Ours 


Fig. 7. When processing an input with significant density variation - e.g., from four synthetic scans (most-left), 
FSSR and ours can avoid generating unwanted artifacts caused by high frequency noises. The total time of our 
reconstruction is 6.81 sec. (s = 3.0 and 342k triangles are obtained on the resultant mesh), while FSSR takes 
670 sec. and results in 301 k triangles (scale=0.0105). Both are tested on a CPU with eight-cores. 


as the scale used in FSSR. This is consistent with the 
formulation presented in [5], where the support size 
is set as three times of the input scale. The actual 
values of scale and sd are also given in Fig.6. It can be 
found that similar number of triangles are generated 
in FSSR and our method by setting the value of scale 
in this way. Note that, as both FSSR and ours do not 
need any global operations during the computation 
of surface reconstruction, it can easily be parallelized 
on the PC with multi-cores - OpenMP is used in our 
implementation. We test both approaches on a PC 
with 8-cores. As shown in the computational statistics 
in Fig.6, the program can be effectively speed up on 


8-cores. The multi-core version of FSSR is provided 
by the authors on their homepage. 

We also study the effectiveness of our approach on 
the benchmark of FSSR with shape density variation 
caused by superposing point sets obtained from mul¬ 
tiple scans. As shown in Fig.7, the input set from four 
synthetic scans is downloaded from the homepage of 
FSSR's authors. Our reconstruction is similar to the 
result from FSSR but ours method is much faster. 

5.3 Verification of numerical error 

Error bound of the quasi-solution A obtained by our 
closed-form formulation with reference to the exact 
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Fig. 8. Examples of reconstruction from sets (having 250k points) with different level of Gaussian noises. 
Topological errors can be found on the results of Screen-Poisson and FSSR (see the regions circled by dashed 
lines in red). The cross-sectional view of function values in our reconstruction has also been given in the right, 
where the regions in white color have undefined function values. 


TABLE 4 

Error Statistics of Quasi-Solution 


Model 

Figure 

V 

||A-AH* 

Ramesses 

5(a) 

457,616 

9.52 x 10~ s 

Raptor 

5(b) 

1,666,700 

1.98 x 10“ 8 

Aquarius 

6(a) 

176,771 

3.46 x 10~ 7 

Horse 

6(a) 

149,459 

3.47 x 10 -7 


solution A of Eq.(7) has been derived in Section 3.2. 
It is also interesting to study the error between A and 
A. We measure ||A — A x in examples shown above 
and the results are listed in Table 4. 

From the statistics, it can be easily found that our 
quasi-solution provides very accurate results on both 
the clean point cloud and the real data. The numerical 
solver for computing the exact solution runs out of 
memory on the two examples - Momenta and Nep¬ 
tune in Fig.5. Thus, the errors cannot be evaluated and 
shown here. 

5.4 Noisy data 

In the following tests, we verify the robustness of our 
approach on input with noises at different levels. For 
a given point set V with normal vectors A f, if the 
diagonal length of its bounding box is d, a new point 
set with 1 5% Gaussian noises is obtained as follows: 

• ng = r j^n'p-' points are selected from V into a 
sub-set Q with n-p denoting the number of points 
in V; 

• Randomly generate a set of scale with Gaussian 
distribution, T>g = {d;}, i = 1,2, ...,nc, with 
di e [0, (M/1000]; 

. Impose the noises onto the points in Q by p' = 
p y + djUj for all p, G Q. 



Fig. 9. Color maps to visualize the forward-distance 
based errors on the results generated by FSSR (top 
row) and our method (bottom row). 

Normal vectors on a set of noisy points are re¬ 
generated by the orientation-aware Principal Compo¬ 
nent Analysis (PCA) with p-nearest neighbors - here 
p = 6 is used in all tests. 

We reconstruct mesh surfaces from a filigree model 
with 30% and 60% Gaussian noises by different meth¬ 
ods, including Screen-Poisson, FSSR and ours (see 
Fig.8). The noise-free set of points are sampled from 
a mesh model so that we can compare the results 
of reconstruction with the original mesh to evaluate 
the shape approximate errors generated by different 
method. Screen-Poisson reconstruction does not per¬ 
form well on noisy model. As shown in Fig.8, models 
with incorrect topology are generated. FSSR and ours 
can still reconstruct 'correct' models even after embed¬ 
ding 60% Gaussian noises. We then compare these two 
methods in terms of shape approximation error by 
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TABLE 5 

Error Statistics of Reconstruction on Noisy Input 


j FSSR on the Filigree Model (Fig.8) | 

Noisy 

Level 

Forward Dist. 

Backward Dist. 

Scale 

Max. 

Ave. 

Max. 

Ave. 

10% 

.00360 

.000295 

.0175 

.00315 

.00375 

30% 

.00800 

.00600 

.0325 

.0014 

.00670 

60% 

.0115 

.00145 

.0300 

.00170 

.0112 

j Our Method on the Filigree Model (Figi 

9 

Noisy 

Forward Dist. 

Backward Dist. 


Level 

Max. 

Ave. 

Max. 

Ave. 

S 

10% 

.00405 

.000255 

.0115 

.000275 

1.9 

30% 

.00700 

.000800 

.00700 

.000800 

2.7 

60% 

.0135 

.00210 

.0150 

.00220 

3.5 


*The errors are measured by the Metro tool [41]. 


the Metro tool [41] (see Table 5). In the measurements 
based on forward distances from ground-truth to the 
reconstruction, FSSR has slightly smaller errors (see 
also the visualization in Fig.9). In the errors based 
on backward distances (i.e., from reconstruction to 
ground-truth), our method outperforms FSSR. This is 
because FSSR generates some interior isolated regions 
(i.e., topological errors) but our method does not - 
see the zoom-view in Fig.8. Moreover, our method is 
17.5x and 36.4x faster than FSSR on the 30% and 60% 
noisy models respectively. 

5.5 Limitation 

The limitations of our approach are mainly caused 
by the nature of locally compact support of kernel 
functions. As a result, we share the following common 
limitations as the FSSR method. 

• Near the boundary of regions with function-value 
defined, some small fragments isolated from the 
main reconstruction could be formed by the nu¬ 
merical oscillation. Such isolated fragments must 
be removed by the post-processing step taken on 
the mesh surface after polygonization. 

• Although reconstruction with high quality can be 
found at the example shown in Fig.7, misaligned 
multiple scan could lead to multi-layers of points, 
therefore also have multiple surface layers pro¬ 
duced at those 'overlapped' regions. 

Caused by these limitations, when a set of points 
with low quality (e.g., the set obtained from two 
Kinect sensors as shown in Fig.10) are used as input, 
neither FSSR nor ours can obtain water-tight surface 
as generated by Screen-Poisson reconstruction. 

6 Conclusion 

In this paper, we present a novel surface reconstruc¬ 
tion method based on computing an approximate 
solution of FIRBF-based implicit surface fitting. The 
approximate solution is formulated as a weighted 
sum of compactly supported basis functions centered 



Fig. 10. For an input point set (left most) with low 
quality from Kinect, all methods generate poor results. 

at input data points equipped with normal vectors 
(i.e., we provide a closed-form solution without any 
global operation). The implicit function for surface re¬ 
construction can be efficiently and robustly evaluated. 
The error-bound between our approximate solution 
and the exact solution has been derived, which can 
be guaranteed as long as the maximum number of 
points covered in the supports of RBF kernels is 
capped by a fixed number. Moreover, to strengthen 
the performance of our approach on input with high 
non-uniformity, a center selection algorithm has also 
been introduced. Experimental results have shown the 
performance of our approach by comparing to the 
state-of-the-arts. 

No global operation needs to be applied during the 
surface reconstruction of our approach. As a result, 
it is easy to extend our implementation to run in the 
out-of-core manner or on a distributed PC-cluster. We 
would like to further investigate the strength of our 
method in this aspect in our future work, which can 
make it possible to realize the on-site reconstruction 
of large-scale 3D models (e.g., outdoor scenes like city 
scale). Many robotics and virtual reality applications 
could benefit from this work. 
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